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Preface 


Interest  in  a  spacecraft's  orbital  trajectories 
about  the  Martian  moon,  Phobos,  has  arisen  with  renewed  U.S. 
and  Soviet  emphasis  on  exploration  of  the  Martian  system. 

The  Report  of  the  National  Commission  on  Space,  appointed  by 
President  Reagan  and  chaired  by  Thomas  Paine,  recommended  a 
thorough,  efficient  and  systematic  progression  towards  Mars 
(11:193)  to  give  NASA  focus  and  revive  a  sagging  U.S.  space 
program.  The  report  recommends  unmanned  probes, 
penetrators,  and  sample  return  missions  to  the  Moon,  to  Mars 
and  its  moons,  to  some  promising  asteroids,  and  to  the  outer 
planets  and  their  moons  followed  by  automated  mining  and 
materials  processing  plants,  and  eventually  by  manned 
explorations  and  human  outposts.  The  Planetary  Society  is 
pushing  for  the  U.S.  to  make  an  official  declaration  to 
strive  towards  human  exploration  of  Mars  (12:3).  The 
Soviets  have  already  stated  they  are  pursuing  the 
possibility  of  landing  a  man  on  Mars  (11:161;  6:14-15). 
Soviet  records  for  long  duration  space  flight  set  in  their 
earth  orbiting  space  station  provide  groundwork  for  a  long 
duration  manned  flight  to  Mars.  The  Soviets  currently  have 
two  spacecraft  on  their  way  to  the  Martian  system  and  have 
another  Mars  mission  scheduled  for  1994  (6:15;  2:16;  4:392). 
The  next  approved  U.S.  mission  to  Mars,  the  Mars  Observer, 
is  scheduled  to  launch  in  1992  (4:392). 
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The  two  Soviet  spacecraft,  Phobos  I  and  Phobos  II, 
were  launched  July  7  and  July  12  of  1988  and  should  be 
arriving  in  the  Martian  system  near  the  end  of  January  1989 
(4:392).  The  480  million  dollar  (10:9b)  mission  includes  a 
lander  carried  on  each  spacecraft  and  Phobos  II  carries  a 
surface  hopper  intended  to  land  and  probe  the  surface  of 
Phobos  (4:392-393:  5:183).  With  some  success,  Phobos  could 
"become  the  fourth  extraterrestrial  body  on  which  spacecraft 
from  earth  have  landed ."( 4 : 392 )  On  September  2,  a  Soviet 
ground  controller  sent  Phobos  I  an  incorrect  command  causing 
loss  of  attitude  control  of  the  spacecraft  and  its  solar 
panels  (5:183).  With  the  solar  panels  improperly  aligned, 
there  was  not  enough  energy  to  sustain  the  transmitter. 
Soviet  officials  are  hoping  the  solar  panels  will  get 
aligned  with  the  sun  and  restore  power.  Early  Soviet 
missions  to  Mars,  all  before  1974,  did  not  meet  with  great 
success  (4:392).  They  either  crashed,  missed  their  target 
or  stopped  transmitting  data  early. 

The  last  earth  vehicles  to  go  to  Mars  were  the 
successful  U.S.  Viking  orbiters  and  landers  which  reached 
+he  planet  in  the  summer  of  1976  (4:392).  Viking  flybys  of 
Phobos  provided  the  estimate  of  its  gravitational  parameter 
used  in  this  study  (13:35). 
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Abstract 

Orbital  trajectories  obtained  from  numerical 
integration  of  over  two  thousand  sets  of  initial  conditions 
for  equations  of  motion  of  a  spacecraft  in  the  Mars  Phobos 
system  are  examined  in  a  Phobos  centered  rotating  coordinate 
frame.  The  equations  of  motion,  simplified  with  the  choice 
of  the  system  of  units,  are  integrated  using  a  Hamming 
fourth  order  predictor  corrector  algorithm.,  The 
trajectories  were  examined  by  plotting  the  position  of  the 
spacecraft  and  by  listing  of  the  state  vector  values  at  each 
crossing  of  the  X,  Y,  and  Z  axes. 

."As  initial  velocity  and  altitude  are  varied, 
trajectories  in  the  orbital  plane  of  Phobos  about  Mars 
follow  an  orderly  pattern.  A  range  of  initial  velocities 
resulted  in  orbital  trajectories  about  Phobos  at  a  given 
altitude.  Near  the  center  of  this  range  of  initial 
velocities,  termed  the  orbital  window,  is  a  unique  initial 
velocity  that  resulted  in  a  closed  periodic  orbit  (at  a 
given  altitude)  in  the  Phobos  centered  rotating  coordinate 
frame.  Initial  velocities  greater  than  or  less  than  the 
velocity  needed  for  a  closed  periodic  orbit  result  in 
trajectories  that  move  away  from  the  periodic  orbit  in  a 
predictable  manner  and  eventually  leave  the  orbit  window. 
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Trajectories  outside  the  orbit  window  either  collide  with 
Phobos  or  leave  the  vicinity  of  Phobos . 

Trajectories  out  of  the  orbital  plane  of  Phobos 
about  Mars  are  more  complex,  three  dimensional  (helical 
shaped)  paths  that  do  not  remain  in  an  inclined  plane,  but 
do  exhibit  some  order.  Again  there  is  a  range  of  initial 
velocities  that  fall  within  an  orbital  window.  No  orbital 
trajectories  were  found  that  remained  in  a  plane  which 
contained  Phobos  and  was  inclined  to  the  orbital  plane  of 
Phobos  about  Mars.  No  closed  periodic  orbits  were  found 
outside  the  orbital  plane  of  Phobos  about  Mars. 


NUMERICAL  STUDY  OF  ORBITAL  TRAJECTORIES  ABOUT  PHOBOS 


I .  Introduction 


Background 

Phobos  is  the  larger  of  two  small  moons  orbiting 
Mars.  Phobos  has  been  described  as  a  potato  shaped  rock 
(4:392)  and  has  been  modeled  mathematically  as  a  triaxial 
ellipsoid  by  Werner  (13:1).  Phobos  is  a  gravity  gradient 
stabilized  satellite  with  its  long  axis,  length  of  27  ±  1  Km 
(13:35),  maintained  in  a  Mars  pointing  orientation  by  the 
gravity  gradient  torques.  Phobos ' s  length  along  the  axis  in 
its  direction  of  motion  is  21.6  ±  1.4  Km  and  its  length 
along  the  axis  out  of  its  orbit  plane  is  18.8  ±  1.4  Km 
(13:35).  The  gravitational  parameters  for  Phobos  and  Mars 
are  6.6e-4  and  42826.32  Km3/sec2  respectively  (13:35). 
Because  the  mass  of  Phobos  is  about  65  million  times  less 
than  the  mass  of  Mars,  Mars  is  the  dominating  force  on  a 
spacecraft  orbiting  in  the  Mars  Phobos  system.  Phobos  is  in 
a  circular  orbit  9378  Km  from  the  center  of  Mars  and 
completes  an  orbit  in  7.65  hours  (13:35). 

Werner  used  Hamiltonian  mechanics  to  develop 
restricted  three  body  equations  of  motion  for  a  spacecraft 
in  orbit  in  the  Mars  Phobos  system  (13:1-12).  He  then 
searched  for  and  found  closed  periodic  orbits  about  Phobos 
in  the  plane  of  Phobos ' s  orbit  about  Mars.  He  showed  these 
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were  unstable  orbits  and  suggested  the  existence  of 
bifurcation  regions  caused  by  nearby  inclined  orbits  but 
found  no  inclined  periodic  orbits  that  did  not  pass  through 
Phobos  (13:19-20).  A  bifurcation  is  "where  a  sudden  change 
in  behavior  occurs  as  a  parameter  passes  through  a  critical 
value"  (8:317).  It  was  conjectured  that  the  equations  of 
motion,  being  nonlinear,  could  give  rise  to  a  variety  of 
interesting  orbital  trajectories  as  some  initial  conditions 
or  control  parameters  are  varied. 

Werner's  numerical  exploration  of  his  equations  of 
motion  was  limited  to  finding  closed  periodic  orbits  in  the 
plane  of  Phobos ' s  orbit  about  Mars  and  discussing  their 
stability  using  Floquet  theory.  Additional  work  was 
necessary  to  classify  a  broader  spectrum  of  the  possible 
orbital  trajectories  about  Phobos  both  in  and  out  of  the 
plane  of  its  orbit  about  Mars. 

Problem  Statement 

The  intent  of  this  research  effort  is  to  begin  with 
Werner's  equations  of  motion,  develop  and  apply  computer 
software  to  explore,  examine,  and  classify  the  possible 
orbital  trajectories  near  Phobos  both  in  and  out  of  its 
orbital  plane  and  to  look  for  possible  bifurcations. 
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Approach 


The  general  approach  was  to  numerically  integrate 
the  equations  of  motion  beginning  with  a  large  number  of 
different  sets  of  initial  conditions  chosen  in  a  way  that 
would  yield  a  set  of  orbital  trajectories  adequately 
representing  all  the  possible  trajectories.  Orbital 
trajectories  in  the  plane  of  Phobos ' s  orbit  about  Mars  were 
obtained  by  integrating  over  1500  different  sets  of  initial 
conditions.  Orbital  trajectories  out  of  the  plane  of 
Phobos 1 s  orbit  about  Mars  were  obtained  by  integrating  over 
700  different  sets  of  initial  conditions. 

The  trajectories  obtained  by  these  numerical 
integrations  were  examined  using  plots  and  listings 
extracted  from  the  state  vector  time  history.  Plots  for  a 
few  hundred  of  the  trajectories  were  produced  from  samples 
of  the  state  vector.  The  state  vector  was  sampled  at 
approximately  44  second  intervals  to  keep  the  size  of  the 
plot  files  small  while  using  enough  data  for  the  plots  to 
appear  continuous.  Listings  produced  for  all  of  the 
trajectories  displayed  the  state  vector  at  each  crossing  of 
the  X,  Y,  or  Z  axes,  the  period  about  Phobos  in  the  XY 
plane,  and  messages  to  inform  when  a  closed  orbit  was 
completed,  when  a  collision  occurred,  or  when  the 
distance  from  Phobos  exceeded  the  Mars-Phobos  distance. 

The  set  of  possible  trajectories  were  separated  into 
those  that  remained  in  Phobos ' s  orbital  plane  about  Mars, 
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and  those  that  did  not.  These  were  further  divided  into 
trajectories  that  circled  Phobos  at  least  once  (orbital)  and 
those  that  did  not. 

Closed  periodic  orbits  about  Phobos  were  special 
cases  of  the  orbital  trajectories  that  were  investigated.  A 
closed  periodic  orbit  was  one  that  later  returned  to  the  set 
of  initial  conditions  for  position  and  velocity.  A  closed 
orbit  was  detected  by  the  software  if  the  state  vector 
returned  to  the  initial  conditions  to  within  a  set  of 
tolerances  after  circling  Phobos  at  least  once.  This  set  of 
tolerances  was  defined  as  the  set  containing  the  largest 
change  in  each  state  variable  encountered  during  any 
integration  step  since  the  beginning  of  the  integration  for 
a  particular  trajectory. 

XY  Planar  Trajectories 

The  initial  values  of  the  Y  coordinate  (altitude) 
and  the  X  velocity  are  varied  to  generate  families  of 
orbital  trajectories  in  the  XY  plane  (the  plane  of  Phobos ' s 
orbit  about  Mars).  All  the  other  initial  conditions  are 
taken  to  be  zero.  With  these  initial  conditions,  the  total 
spacecraft  velocity  is  simply  the  initial  value  of  the  X 
velocity.  Therefore,  the  initial  value  of  the  X  velocity 
for  a  given  initial  +Y  coordinate  is  chosen  as  the  control 
parameter  to  generate  orbital  trajectories  about  Phobos  in 
the  XY  plane.  Zero  values  for  the  initial  Z  coordinate  and 
the  initial  Z  velocity  restrict  the  motion  of  the  satellite 
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to  the  XY  plane  since  each  term  in  the  equations  for  Z  and  Z 
contains  a  Z  or  a  Z.  The  choice  of  zero  initial  values  for 
the  X  coordinate  and  the  Y  velocity  makes  the  initial 
velocity  of  the  satellite  simply  the  chosen  control 
parameter,  the  X  velocity,  which  is  then  tangent  to  the 
trajectory.  Choosing  the  initial  conditions  this  way 
enabled  thinking  of  the  selection  of  an  orbital  trajectory 
as  a  choice  of  the  kinetic  energy  (determined  by  the  choice 
of  the  X  velocity)  and  the  potential  energy  (determined  by 
the  choice  of  the  Y  coordinate).  Then,  the  problem  of 
selecting  a  particular  type  of  orbit  about  Phobos  could  be 
related  to  the  simpler  two-body  problem  of  orbit  selection. 
Given  a  particular  altitude  in  the  two-body  problem,  a 
certain  velocity  generally  defines  a  closed  orbit.  However, 
in  the  restricted  three-body  problem,  given  a  Y  coordinate 
(altitude),  selection  of  the  X  velocity  defines  an  orbital 
trajectory  which  in  general,  is  not  a  closed  path. 

Three  Dimensional  Trajectories 

To  generate  trajectories  out  of  the  XY  plane,  the 
initial  values  of  the  Z  coordinate  and  the  X  velocity  were 
varied  for  a  given  initial  Y  coordinate  value.  Again,  since 
the  initial  values  of  the  Y  velocity  and  the  Z  velocity  were 
zero,  the  total  spacecraft  initial  velocity  was  simply  the 
initial  X  velocity  which  again  could  be  treated  as  the 
control  parameter.  In  an  attempt  to  find  indications  of  a 
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bifurcation,  the  initial  Y  coordinate  of  20  Km  was  chosen 
because  it  fell  within  a  region  described  as  a  possible 
bifurcation  region  by  Werner  (13:19). 

Overview 

The  introduction  presented  in  this  chapter  gives 
some  background  information  leading  to  the  problem 
statement,  then  gives  the  general  approach  taken  to  solve 
the  problem. 

Chapter  II  begins  the  problem  development  with 
Werner’s  restricted  three  body  equations  of  motion  for  a 
spacecraft  in  the  Mars  Phobos  system.  The  Phobos  centered 
rotating  coordinate  frame  of  reference  used  for  these 
equations  is  described.  A  system  of  units  for  mass, 
length,  and  time  is  chosen  which  simplify  the  equations  of 
motion . 

Chapter  III  describes  the  method  of  solution.  To 
solve  the  problem  of  describing  the  possible  orbital 
trajectories  about  Phobos,  software  was  needed  to  enable 
examination  of  a  large  number  of  numerical  solutions  to 
equations  of  motion.  The  software  developed  to  integrate 
the  equations  of  motion  from  many  sets  of  initial  conditions 
for  the  spacecraft's  position  and  velocity  is  provided  in 
the  appendix. 

The  results  are  presented  in  Chapter  IV  using  graphs 
that  divide  the  solution  space  into  several  categories 
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representing  different  outcomes  for  the  trajectories 
obtained  from  the  integrations.  A  number  of  plots  are 
presented  showing  typical  trajectories  that  fall  into  these 
categories.  The  results  are  discussed  in  three  sections 
that  separate  the  trajectories  into  those  that  result  in 
closed  periodic  orbits  in  the  XY  plane,  XY  planar 
trajectories  in  general,  and  three  dimensional  trajectories. 

Chapter  V  presents  conclusions  about  practical 
implications  when  orbiting  a  spacecraft  about  Phobos  which 
are  drawn  from  the  results  of  Chapter  IV.  They  include 
discussion  of  the  amount  of  precision  in  control  of  the 
state  vector  needed  and  the  amount  and  frequency  of  velocity 
adjustments  needed  to  stay  in  the  orbit  window. 
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I I .  Problem  Development 


The  motion  of  a  spacecraft  of  negligible  mass 
orbiting  in  the  Mars  Phobos  system,  a  restricted  three  body 
problem,  was  described  mathematically  by  Werner  (13:11-12). 
He  modeled  the  attraction  due  to  the  irregular  shape  of 
Phobos  using  the  potential  energy  field  of  a  homogeneous 
triaxial  ellipsoid  rather  than  a  sphere  which  is  usually 
treated  as  a  point  source.  Using  Hamiltonian  mechanics,  he 
derived  the  following  equations  of  motion  for  the  system. 


X  =  P^  +  Q(Y  -  D) 


(1) 


Y  *  Py  -  QX 


(2) 


Z  =  P, 


(3) 


P„  =  QPy  -  GM0(X/D3  -  3XY/D*  +  3RaX/2D°)  -  GMaX/R3 

+  3GIX/4R°  -  3GX/4R7[(5Xa  -  2Ra)(I-2Iw) 


+  5Y  (I  -  Iw)+  5Z  ( I  -  2I„)  ]  (4) 


Pv  =  -  QP„  +  GM0[1/Da  -  4Y/D3  -  3 ( YRa  -  DRa  - 

2DYa)/2D8]  -  GMjY/R3  +  3GIY/4R®  - 
3GY  [  ( 5Ya  -  2Ra )  ( I  -  2IW)  +  5Xa(I  -  2IW)  + 

5Za  ( I  -  2I«)]/4R7  (5) 
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P«  =  -  GM0(Z/D3  -  3YZ/D*  +  3ZRa/2D°)  -  GMjZ/R3  + 

3GZI/4R°  -  3GZ[(5Za  -  2Ra)(I  -  21.*)  + 

5X2(I  -  2Ira)  +  5Ya(I  -  2Iyy)]/4R?  (6) 

Where  X,  Y,  Z  are  the  coordinates  of  the  spacecraft  position 
(altitude)  and  X,  Y,  Z  are  the  components  of  the  spacecraft 
velocity  expressed  in  the  coordinate  frame  described  below. 
P„,  Py<  and  P„  are  the  conjugate  momenta  (9:172-173). 

The  other  terms  used  in  eqs  ( 1 ) — ( 6 )  are 

P  =  angular  velocity  of  XYZ  coordinate  frame 
Ra  =  Xa  +  Ya  +  Za 
M0  =  mass  of  Mars 

Mi  =  mass  of  Phobos 

G  =  universal  gravitational  constant 

D  =  Mars  to  Phobos  distance 

I  =  +  Iw  +  I„  =  sum  of  the  principal 

moments  of  inertia  of  Phobos 
The  principal  moments  of  inertia  of  a  homogeneous 
ellipsoid  are  (7:501) 

I«c  =  Mi  (ba  +  ca)/5  (7) 

Iw  =  Mi  ( a*  +  ca)/5  (8) 

I«  =  Mj.  ( aa  +  ba)/5  (9) 

Where 

a  =  half  length  of  Phobos  along  X  axis 

b  =  half  length  of  Phobos  along  Y  axis 

c  =  half  length  of  Phobos  along  Z  axis 
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The  coordinate  system  for  eqs  ( 1 ) — ( 6 )  is  a  rotating 
cartesian  system  with  the  positive  Y  axis  always  pointed 
towards  Mars  and  the  positive  X  axis  pointed  in  the 
direction  of  Phobos  motion  in  its  orbit  about  Mars.  Figure 
1  shows  Mars  and  Phobos  in  the  XY  plane.  The  size  of  Phobos 
is  scaled  up  10  times  to  make  it  just  visible.  The  Z  axis 
points  out  of  the  plane  of  the  orbit  of  Phobos  (out  of  the 
page)  forming  a  right-handed  coordinate  system. 
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Figure  1.  Phobos  Centered  Rotating  Coordinate  System 
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The  choice  of  the  system  of  units  generally  used 
with  restricted  three  body  problems  allows  some 
simplification  of  the  equations  of  motion.  The  mass  unit  is 
defined  by  letting  the  sum  of  the  masses  of  the  two  primary 
bodies  be  equal  to  1  mass  unit  (M.U.).  The  unit  of  length, 

1  distance  unit  (D.U.),  is  the  distance  between  the  two 
primaries  (9378  Km).  The  time  unit  (T.U.)  is  chosen  so  that 
the  gravitational  constant  G  is  equal  to  1.  This  happens 
when  the  period  of  the  two  primaries  in  their  orbit  about 
each  other  (7.65  hours)  is  equal  to  2tt  T.U.  One  T.U.  is 
then  equal  to  4383  seconds  (about  1.2  hours).  The  angular 
velocity,  Q,  of  the  coordinate  frame  is  then  1  radian/T.U. 
With  these  units,  and  letting 

1 1  ~ I»oc  "t  1  yy  "t  Iu 

la  “  i  >oc  —  lyy  + 

I.  »  I«  +  Iw  -  I  — 

where  the  principal  moments  of  inertia  are  calculated  using 
eqs  (7) -(9)  with  the  just  defined  units  for  the  mass  and 
length  terms,  Eqs  (1)  -  (6)  simplify  to 


X 


P*  +  Y  -  1 


(10) 


Y  =  Py  -  X  (11) 

Z  =  P«  (12) 

Px  =  Py  -  M0(X  -  3XY  +  3R*X/2 )  -  MJC/R3  +  3IX/4R®  - 

3X[(5X“-  2Ra) I i  +  5YaIa  +  5ZaI3]/4R7  (13) 
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-  P~  +  Mc[l  -  4Y  -  3 ( YRa  -  Ra  -  2Ya)/2]  -  MjY/R3 

+  3IY/4R°  -  3Y[(5Y2  -  2Ra)Ia  +  5XaIx  +  5ZaI3]/4R7  (14) 

-  M0(Z  -  3YZ  +  3ZRa/2 )  -  MjZ/R3  +  3ZI/4R®  - 

3Z[(5Za  -  2Ra ) I 3  +  5X*!,.  +  5Y2I2]/4R7  (15) 
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Ill . 


Method  of  Solution 


Eqs  (10)— (15)  are  the  coupled  nonlinear  differential 
equations  to  integrate  to  get  the  path  of  a  spacecraft  in 
the  Mars  Phobos  system.  The  integration  is  performed  by  the 
FORTRAN  program  listed  in  the  Appendix.  Eqs  (10)— (15) 
appear  in  ^he  subroutine  "RHS.for",  called  by  subroutine 
"Haming . f or " ,  a  fourth  order  predictor  corrector  integration 
algorithm.  The  main  program,  "Phobos . for " ,  reads  a  set  of 
initial  conditions  and  the  length  of  time  to  integrate  from 
an  input  file,  "IC.dat".  The  program  integrates  the 
equations  of  motion  until  the  trajectory  terminates  in  a 
collision  with  Phobos  or  Mars,  exceeds  1  D.U.  from  Phobos, 
or  the  specified  length  of  time  is  reached.  Then  a  new  set 
of  initial  conditions  are  read  from  the  input  file  and  the 
program  integrates  again  and  again  until  the  end  of  file  is 
read  from  "IC.dat". 

The  program  types  selected  information  to  the  screen 
(which  can  be  saved  in  a  log  file)  and  writes  the  state 
vector  (position  and  velocity)  and  the  elapsed  time  to  data 
files  if  plots  are  desired.  The  information  typed  to  the 
screen  includes  a  notice  when  a  coordinate  is  passing 
through  a  value  of  zero,  the  rate  of  change  of  that 
coordinate,  and  the  value  of  the  X  coordinate  (or  the  Y 
coordinate  if  X  is  passing  through  zero).  Also  typed  to  the 
screen  are  the  time  to  complete  each  orbit,  the  orbit  count, 


and  a  message  to  indicate  the  reason  for  early  termination 
of  the  integration  if  needed. 

The  two  options  to  write  the  state  vector  and 
elapsed  time  to  data  files  allows  the  data  to  be  plotted 
using  any  plotting  software  that  plots  columns  of  data 
against  another  column  (like  a  spreadsheet).  The  first 
option,  to  write  the  state  vector  and  the  elapsed  time  to  an 
output  file  named  " state . dat ; n"  each  .01  T.U.  (43.83  sec), 
enabled  the  plotting  of  the  orbital  trajectories.  The 
second  option,  to  write  the  same  information  to  a  file  named 
"sect  ion . dat ; n"  each  time  the  trajectory  crosses  the  Y  =  0 
plane,  is  intended  to  produce  plane  section  plots. 

The  integration  step  size  was  chosen  large  enough  so 
that  the  computer  run  time  for  a  group  of  trajectories  was 
reasonably  short,  yet  small  enough  so  that  a  smaller  step 
size  wouldn't  significantly  change  the  output  state  vector 
or  the  accumulated  period  of  an  orbital  trajectory.  A 
significant  change  in  the  state  vector  was  one  that  was 
larger  than  the  largest  change  encountered  during  any 
integration  step  for  the  orbit.  A  significant  change  in  the 
orbital  period  was  one  that  was  larger  than  the  step  size 
chosen.  The  choice  of  0.0001  T.U.  (.4383  seconds)  for  the 
step  size  accomplished  this  goal.  A  10  T.U.  trajectory 
would  run  in  about  a  minute.  Using  a  0.00001  T.U.  step  size 
didn't  change  the  results  but  took  much  more  computer  time. 

A  closed  orbit  is  detected  by  the  software  when  the 
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state  vector  returns  to  the  initial  conditions  within  an 
amount  equal  to  the  largest  step  size  of  the  state  vector 
during  the  integration.  Because  that  step  size  depends  on 
the  integration  time  step,  a  small  time  step  is  desirable. 
The  time  reported  for  the  orbital  period  also  depends  on  the 
step  size.  The  chosen  step  size,  less  than  0.5  seconds . 
makes  the  reporting  of  the  period  to  the  nearest  second 
simpler  and  more  convenient  than  a  larger  step  size  which 
would  have  required  some  interpolation  between  steps  to 
report  to  that  accuracy. 

The  trajectories  are  not  allowed  to  pass  through  the 
surface  of  Phobos  or  Mars.  The  position  of  the  spacecraft 
is  checked  each  step  of  the  integration  to  ensure  it  falls 
outside  the  equation  of  an  ellipsoid  centered  on  the  origin 
with  Phobos ' s  dimensions.  The  integration  is  stopped  if  the 
position  of  the  spacecraft  falls  inside  the  ellipsoid  for 
three  reasons.  First,  the  equations  of  motion  are  no  longer 
valid  inside  the  surface  of  Phobos  because  the  potential 
energy  expression  used  to  develop  the  equations  is  different 
inside  the  ellipsoid  (3:08-107).  Second,  numerical 
difficulties  can  arise  if  the  position  of  the  spacecraft  is 
rear  the  origin  (R  approaching  zero)  because  the  equations 
are  full  of  terms  divided  by  P.  Third,  the  spacecraft 
crashes  into  Phobos  and  the  flight  is  over.  The  position  of 
vhe  spacecraft  is  also  checked  to  ensure  it  remains  outside 
a  spher*-  with  a  Mars  radius  centered  on  V  equals  1  D  .  ”  . 
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because  that  also  represents  an  undesirable  trajectory  that 
crashes  into  Mars. 

No  trajectories  were  encountered  that  ended  in 
collision  with  Mars  because  all  the  trajectories  considered 
were  near  Phobos  and  therefore  had  plenty  of  orbital 
velocity  about  Mars. 
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IV.  Results 


CIos ed  Periodic  Crbi ts  In  the  XY  Plane 

Werner  found  closed  periodic  orbits  are  possible  in 
the  XY  plane  from  near  the  surface  of  Phobos  to  beyond  5000 
Km  (13:23).  With  proper  selection  of  the  initial  X  velocity 
(the  chosen  control  parameter),  a  closed  orbit  in  the  XY 
plane  can  be  obtained  for  any  practical  altitude.  All  these 
closed  planar  orbits  proceed  in  a  clockwise  direction  as 
seen  looking  down  the  +Z  axis.  Figure  2  shows  some  closed 
periodic  orbits  about  Phobos  with  Y  axis  crossings  at  20, 

40,  60,  80,  100,  125  Km  altitude. 

Figure  3  shows  that  for  periodic  orbits  the  value  of 
the  X  velocity  is  generally  positive  (spacecraft  velocity 
greater  than  Phobos  velocity)  for  Y  greater  than  zero 
(closer  to  Mars)  and  is  always  negative  (spacecraft  velocity 
less  than  Phobos  velocity)  for  Y  less  than  zero  (farther 
from  Mars).  This  is  the  expected  result  for  a  spacecraft  in 
a  purely  two-body  orbit  about  Mars.  The  velocity  of  the 
orbit  decreases  with  increasing  distance  from  Mars  as 
kinetic  energy  is  traded  for  potential  energy.  The  values 
of  the  X  velocity  for  orbits  with  higher  altitudes  about 
Phobos  decrease  some  because  they  are  in  lower  energy 
elliptical  orbits  about  Mars  than  the  circular  orbit  of 
Phobos  about  Mars.  The  result  is  the  bowed  look  of  the  plot 
in  figure  3.  Orbits  in  the  opposite  direction 
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Xdot  CD 


Figure  3.  X  Velocity  for  Periodic  Orbits  in  the  XY  Plane 
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( counterc lockwise )  are  not  found  and  are  unlikely  because 
Mars  gravity  and  the  resulting  orbital  motion  about  Mars 
dominates  even  when  fairly  close  to  Phobos .  Figure  4  shows 
for  periodic  orbits,  the  Y  velocity  must  be  greater  than 
zero  (moving  closer  to  Mars)  when  X  is  negative  (spacecraft 
behind  Phobos  in  Mars  orbit)  and  the  Y  velocity  must  be  less 
than  zero  (moving  away  from  Mars)  when  X  is  positive 
(spacecraft  ahead  of  Phobos  in  Mars  orbit). 

Using  figures  3  and  4,  given  any  set  of  X  and  Y 
coordinates,  the  components  of  the  spacecraft's  velocity 
needed  in  the  X  and  Y  directions  to  achieve  a  closed 
periodic  orbit  can  be  approximately  determined.  The  choice 
of  the  particular  X  velocity  along  the  width  of  the  bow  in 
figure  3  for  a  given  Y  value  depends  on  the  altitude  of  the 
orbit  about  Phobos. 

With  the  initial  Y  velocity  set  to  zero,  the  initial 
X  velocity  needed  to  produce  a  closed  periodic  orbit  at  a 
desired  Y  altitude  is  plotted  in  figures  5-7.  Figure  5, 
which  shows  the  initial  values  of  the  X  velocity  for  orbits 
near  Phobos,  reveals  a  noticeable  bend  in  the  curve  around  a 
Y  altitude  of  about  20  Km.  The  slow  bend  in  the  line  of 
periodic  orbits  around  the  20  Km  Altitude  may  be  related  to 
the  appearance  of  the  additional  non-zero  pair  of  Poincare 
exponents  found  by  Werner  at  altitudes  near  20  Km  (13:23). 
Figure  6  extends  the  data  out  to  900  Km  altitude  and  appears 
almost  linear.  Figure  7  extends  the  data  further  to  5500  Km 
and  clearly  shows  the  nonlinearity  of  the  data. 
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Initial  Velocities  (Xdot)  for  Periodic  Orbits 
near  Phobos  in  the  XY  Plane 
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Figure  6 


Initial  Velocities  (Xdot)  for  Periodic  Orbits  out 
to  900  Km  in  the  XY  Plane 
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^er iodic  Orbit  Velocities  about  Phobos 


(TnoiBanci5} 

initial  Velocity  (  xoot  )  in  nvs 


Figure  7.  Initial  Velocities  (Xdot)  for  Periodic  Orbits  out 
to  5500  Km  in  the  XY  Plane 
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The  initial  X  velocity,  X^,,  needed  to  generate  a 
periodic  orbit  at  a  given  Y  altitude  near  Phobos  varies 


nearly  linearly  and  can  be  approximated  by  a  linear 
expression  given  by 

Xp(m/s)  =  . 44  Y ( Km )  +  1.2 

After  a  second  guess  for  the  correct  initial  X  velocity  for 
a  periodic  orbit,  a  few  (one  or  two)  simple  linear 
interpolations  quickly  converges  on  the  correct  value.  The 
interpolation  is  accomplished  by  adjusting  the  initial  X 
velocity  to  eliminate  the  difference  in  the  first  two  X  axis 
crossings  for  the  trajectories  obtained  from  the  integration 
of  the  previous  two  guesses. 

Although  time  is  not  explicit  in  the  equations  of 
motion,  the  orbital  period  of  a  spacecraft  in  a  closed 
periodic  orbit  about  Phobos  is  of  interest.  The  orbital 
period  is  obtained  after  the  integration  of  the  complete 
orbit  path  and  is  the  result  of  summing  the  individual 
integration  time  steps.  The  orbital  period  of  orbits  about 
Phobos  increases  rapidly  from  7310  seconds  (about  2  hours) 
for  orbits  at  15  Km  altitude  to  10240  seconds  (almost  3 
hours)  at  50  Km  altitude.  Figure  8  shows  a  plot  of  the 
orbital  periods  for  closed  periodic  orbits  in  the  XY  plane 
from  15  Km  to  200  Km  and  shows  a  drastic  change  in  slope 
around  50  Km.  Above  50  Km  altitude  the  periods  of  the 
orbits  are  about  3  hours  and  increase  slowly  with  increasing 


a  1 1 i t  ude . 


Figure  9  extends  the  plot  of  orbital  periods  out 


Figure  8 


Orbital  Periods  of  Closed  Periodic  Orbits  Near 
Phobos  in  the  XY  Plane 
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to  5500  Km  and  shows  the  slowly  increasing  trend.  The  sharp 
bend  (big  decrease  in  slope)  around  55  to  100  Km  or  10300 
second  period  is  where  Werner  found  an  additional  non-zero 
pair  of  Poincare  exponents  (13:23). 

Werner  found  the  initial  Y  coordinates  (altitudes) 
needed  to  produce  closed  periodic  orbits  in  the  XY  plane 
with  orbital  periods  in  increments  of  100  seconds  beginning 
with  8200  seconds  and  ending  with  11500  seconds  (13:23). 

The  large  jump  in  altitude  he  shows  from  105  Km  with  a  10400 
second  period  to  1284  Km  with  a  10500  second  period  shows 
the  wide  range  of  altitudes  that  have  approximately  the  same 
orbital  period  as  seen  in  figures  8  and  9. 

Figure  10  shows  how  the  orbital  periods  of  the 
closed  periodic  orbits  in  the  XY  plane  found  by  this  study 
compares  to  those  found  by  Werner.  Because  Werner's  data 
are  given  by  incremental  period  and  the  data  in  this  study 
are  given  by  incremental  altitude,  Werner's  data  are 
adjusted  by  interpolation  to  give  periods  at  incremental 
altitudes.  There  is  some  disagreement  (13%  or  less)  in  the 
data  for  orbits  close  to  Phobos.  This  may  be  due  to 
differences  in  the  details  of  the  numerical  integration  such 
as  the  time  step  used  or  due  to  numerical  round-off 
differences  arising  from  the  different  magnitudes  of  the 
numbers  resulting  from  different  choices  for  units  of 
measure.  The  agreement  is  very  good  for  orbits  beyond  50  Km 
alt i tude . 
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Figure  9.  Orbital  Periods  of  Closed  Periodic  Orbits  c.oout 
Phobos 
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Orbital  Periods  about  Phobos 
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Figure  10.  Orbital  Periods  Compared  with  Werner's 
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XY  Planar  Trajectories 

Figure  11  is  a  map  of  the  type  of  trajectories 
obtained  by  varying  the  initial  X  velocity  (kinetic  energy) 
for  a  given  Y  altitude  (potential  energy) .  The  orbital 
trajectories  that  do  not  form  a  closed  orbital  path  about 
Phobos  eventually  either  collide  with  Phobos  or  leave  the 
vicinity  of  Phobos  (escape).  In  either  case,  the  spacecraft 
can  circle  Phobos  not  at  all,  once,  twice,  or  more  times 
before  collision  or  escape.  The  number  of  orbits  about 
Phobos  before  collision  or  escape  is  related  to  how  close 
the  initial  value  of  the  X  velocity  is  to  the  value  needed 
for  the  periodic  orbit  at  that  altitude. 

Figure  11  shows  there  are  three  regions  of  escape 
velocities.  A  spacecraft  beginning  at  a  Y  coordinate  with 
an  initial  X  velocity  that  falls  in  any  one  of  these  three 
regions  will  leave  the  vicinity  of  Phobos  without  circling 
Phobos . 

Figure  11  also  shows  three  regions  of  velocities 
that  result  in  collisions  with  Phobos  before  completing  an 
orbit  about  Phobos.  One  of  these  regions  is  much  narrower 
than  the  other  two.  It  appears  to  be  a  sliver  split  from  a 
larger  collision  region  by  an  escape  region  which  appears 
only  at  altitudes  above  40  Km. 

The  remainder  of  the  area  on  Figure  11  represents 
trajectories  that  circle  Phobos  at  least  once  and  is  called 
the  orbit  window.  The  line  of  velocities  for  closed 
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Figure  11.  Trajectory  Map  in  Phobos  Orbit  Plane 
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periodic  orbits  of  Phobos  is  at  the  center  of  the  orbit 
window.  The  orbit  window  is  narrow  at  low  Y  altitudes  but 
expands  as  the  altitude  increases. 

Two  notable  rapid  expansions  of  the  orbit  window 
occur  that  stand  out  from  the  general  expansion  trend.  The 
first  is  between  Y  altitudes  of  19  and  25  Km.  The  slope  of 
the  left-hand  boundary  of  the  orbit  window  is  greater  than 
90'“’  between  19  and  20  Km.  The  slope  of  this  same  boundary 
is  less  than  90°  between  23  and  24  Km.  Assuming  the 
boundary  is  continuous,  the  slope  must  be  90°  (infinite)  at 
some  Y  altitude  between  20  and  23  Km.  The  second  rapid 
expansion  of  the  orbit  window  is  between  50  and  60  Km.  The 
slope  of  the  left-hand  boundary  of  the  orbit  window  in  this 
range  appears  to  approach  infinity  also.  These  two  ranges 
of  altitudes  correlate  with  the  two  regions  where  Werner 
found  the  additional  non-zero  Poincare  exponents  and 
suggested  they  indicated  the  existence  of  bifurcation 
regions  in  the  solution  space  of  the  equations  of  motion 
(13:19-20).  They  certainly  signaled  the  existence  of  these 
two  ranges  of  infinite  slopes. 

Figure  12  shows  orbital  trajectories  beginning  at  a 
Y  of  20  Km  altitude  with  initial  X  values  of  -4,  -2,  0, 

2,  4,  6,  8,  and  10  m/s  which  are  all  less  than  the  velocity 
needed  for  a  periodic  orbit.  The  trajectory  with  an  initial 
X  of  -  4  m/s  escapes  and  falls  in  the  first  escape  region 
shown  in  Figure  11.  The  rest  of  the  trajectories  in  Figure 
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Figure  12.  Initial  Velocities  Less  Than  Closed  Periodic 
Orbit  Velocity 


12  all  collide  with  Phobos  without  circling  the  moon.  These 
trajectories  all  fall  in  the  first  collision  region  shown  in 
Figure  1 1 . 

Figure  13  shows  trajectories  beginning  at  a  Y  of  20 
Km  altitude  with  velocities  of  14,  16,  18,  20,  and  22  m/s 

which  are  all  greater  than  the  velocity  for  a  periodic 
orbit.  The  trajectory  that  begins  with  an  initial  velocity 
of  22  m/s  just  barely  misses  Phobos  then  escapes.  This 
trajectory  is  just  over  the  boundary  into  the  third  escape 
region  identified  in  figure  11.  The  rest  of  the 
trajectories  in  figure  13  are  in  the  third  collision  region 
(counting  from  left  to  right)  shown  in  Figure  11. 

Values  of  the  initial  X  velocity  greater  than  Xp 
resulted  in  orbits  whose  X  axis  crossing,  Xc,  shifted  right 


33 


2  0  i  Km 

-  ~'\  /%■. 

/  i  \  /  \\Y\\  / 

y  Phpbos'  '  '  '  x  '  ' 


UUI  -!  1  V 


— f - r  '  -ynT- 

I  i  \  \  \  v 

y  ,  ;  \  v 

!  i  •  ‘  \ !  \ 


'  i  '' — /  '  j  '  II  \ 

\\ ; : :  :14  ;  » \ 

V'l  \  i  ^  16  /'Ttl/0  I  \  , 

T  /  /  i  \i 

1  $  nri/fs  (\ 

1  -v'  /'  >  /  ^ 


"  2-0  irv/s 


-22  m/s 


X  Axis  ( D . U . ) 


34 


.initial  velocity  needed  for  a  periodic  orbit  at  that 

•i  1 1  i  ‘  c  de  .  The  trajectory  ends  in  a  collision  with  Phobos 

after  four  orbits. 

Figure  15  shows  an  orbital  trajectory  about  Phobos 
that  begins  at  a  Y  altitude  of  90  Km  with  an  initial  X 
velocity  of  44  ms.  This  velocity  is  a  little  more  than  the 
initial  X  velocity  needed  for  a  periodic  orbit  at  that 
altitude.  This  time,  the  spacecraft  escapes  after  four 
orbits.  After  the  trajectory  stepped  over  Phobos,  the 
increase  in  the  X  axis  crossings  continued  but  by  decreasing 
amounts  until  the  trajectory  appears  to  orbit  a  point  ahead 
of  Phobos  in  its  orbit  about  Mars. 

Values  of  the  initial  X  velocity  less  than  X^, 
resulted  in  orbits  whose  X  axis  crossing,  Xc,  shifted  left 
(decreased)  by  increasing  amounts  until  the  orbit  either 
collided  with  Phobos  or  escaped.  If  the  change  in  Xc  was 
larger  than  the  X  diameter  of  Phobos,  the  trajectory  could 
step  over  Phobos  and  escape  after  completing  one  or  more 
orbits  of  Phobos.  Then,  the  shift  left  in  the  x  axis 
crossings  continues  but  by  decreasing  amounts  until  the 
trajectory  appears  to  orbit  a  point  behind  Phobos  in  its 
orbit  about  Mars. 


35 


Figure  14.  Velocity  Slightly  More  Than  Needed  for  a 
Periodic  Orbit  Results  in  Collision 


Figure  15.  Initial  Velocity  a  Little  More  Than  Needed  for 
a  Closed  Periodic  Orbit  Resulting  in  Escape 


Three  Dimensional  Trajectories 


Figure  16  is  a  map  of  the  type  of  trajectories 
obtained  by  varying  the  initial  X  velocity  (kinetic  energy) 
and  the  initial  Z  altitudes  with  a  fixed  initial  Y  altitude 
of  20  Km.  The  line  at  Z  equals  zero  represent  the 
trajectories  in  the  XY  plane  at  a  Y  altitude  of  20  Km  as 
shown  in  the  XY  plane  trajectory  map  in  Figure  11. 
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The  trajectory  map  shown  in  Figure  16  is  symmetric 
about  the  line  Z  equals  zero.  There  are  two  regions  where 
the  trajectories  end  in  collision  with  Phobos  before 
>—  o  cir*i  o  r  u  ^  ^  .  These  two  collision  regions  give  a 

third  dimension  to  the  collision  regions  shown  in  Figure  11 
and  discussed  with  the  XY  planar  trajectories.  The  two 
regions  are  ellipsoid  shapes  with  the  larger  ellipsoid 
corresponding  to  the  smaller  X  velocities.  The  other 
regions,  two  escape  regions  and  the  orbit  window, 
represented  at  a  Y  altitude  of  20  Km  in  the  XY  planar 
trajectory  map  of  Figure  11,  are  also  represented  in  Figure 
16  and  thus  given  a  third  dimension. 

The  closed  periodic  orbit  at  the  center  of  the  three 
dimensional  orbit  window  lies  on  a  line  in  the  XY  plane 
trajectory  map  of  Figure  11  but  appears  to  be  a  singular 
point  in  the  trajectory  map  of  Figure  16.  No  closed 
periodic  orbits  are  found  other  than  the  one  at  Z  equals 
zero  . 

Figure  16  shows  a  region  of  trajectories  that  orbit 
once  before  escaping.  These  trajectories  appear  at  Z 
altitudes  greater  than  30  Km.  This  region  is  connected  to  a 
similar  region  that  appears  in  the  XY  plane  trajectory  map 
of  Figure  11  at  Y  altitudes  above  30  Km. 

Figure  16  also  shows  a  narrow  collision  region 
appears  at  Z  altitudes  above  85  Km  with  an  initial  X 
velocity  of  5  m/s. 
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Figure  17  shows  an  orbital  trajectory  projected  or 
the  XY  plane  that  begins  with  an  initial  Y  of  20  Km,  an 
initial  Z  altitude  of  10  Km,  and  an  initial  X  velocity  of  12 
ms.  The  trajectory  ends  in  a  collision  with  Phobos  after 
five  orbits.  Figure  18  shows  the  Z  amplitude  versus  time 
for  this  trajectory  increases  as  does  the  period  of  the 
osci 1 1  at  ion . 

Figure  19  shows  an  orbital  trajectory  projected  on 
the  XY  plane  that  begins  with  an  initial  Y  of  20  Km,  an 
initial  Z  altitude  of  60  Km,  and  an  initial  X  velocity  of  8 
ms.  The  trajectory  ends  in  a  collision  with  Phobos  after 
seven  orbits.  Figure  20  shows  the  Z  altitude  versus  time 
for  this  trajectory  has  a  period  of  about  2rt  T.U.  which  is 
characteristic  of  its  orbit  about  Mars. 
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Cone lus ions 

Mars  gravitational  potential  is  by  far  the  dominant 
force  acting  on  a  spacecraft  orbiting  in  the  Mars  Phobos 
system  even  for  orbits  close  to  Phobos.  It  can  be  said  that 
orbits  about  Phobos  are  merely  orbits  about  Mars  that  happen 
to  go  around  Phobos  with  some  perturbation  effects  due  to 
the  gravitational  potential  of  Phobos. 

If  a  spacecraft  is  to  orbit  Phobos,  it  must  control 
its  velocity  to  stay  inside  the  orbit  window.  The  window  is 
small  at  low  altitudes,  only  a  few  meters  per  second  wide, 
but  widens  as  the  altitude  increases. 

Closed  periodic  orbits  are  available  near  the  center 
of  the  orbit  window  in  the  plane  of  Phobos ' s  orbit  about 
Mars.  A  spacecraft  whose  velocity  can  be  controlled  to 
within  a  few  cm/s  can  achieve  these  closed  periodic  orbits. 
Technically,  these  orbits  are  unstable  because  small  errors 
in  the  state  vector  cause  the  spacecraft  to  drift  away  from 
the  closed  periodic  solution.  But  the  drift  rate  is  slow 
for  orbits  near  the  center  of  the  orbit  window.  Because 
this  drift  rate  is  slow,  the  closed  periodic  solution  is  a 
solution  that  can  represent  a  practical  stability  similar  to 
the  practical  stability  exhibited  by  a  bullet  spinning  about 
its  least  moment  of  inertia  axis.  Although  the  bullet  is 
technically  unstable  spinning  about  that  axis,  the  time  of 
flight  is  small  enough  compared  to  the  amplitude  of  the 
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instability  that  it  is  insignificant.  For  a  spacecraft  in 
orbit  about  Phobos,  small  velocity  adjustments  (a  few  meters 
per  second)  provided  at  relatively  infrequent  intervals,  12 
hours  or  longer  depending  on  altitude,  can  keep  the 
spacecraft  in  the  orbit  window. 

No  closed  periodic  orbit  solutions  were  found 
outside  of  the  plane  of  Phobos ' s  Orbit  about  Mars. 
Nevertheless,  a  spacecraft  can  be  maintained  inside  the 
three  dimensional  orbit  window  with  relatively  small  and 
infrequent  velocity  adjustments  by  picking  regions  away  from 
the  edges  of  the  window. 

The  solution  space  of  trajectories  about  Phobos  is 
well  behaved  and  predictable.  The  search  for  bifurcations 
and  unexplained  chaotic  behavior  failed.  The  two  regions 
where  Werner  suggested  the  appearance  of  additional  Poincare 
exponents  indicated  possible  bifurcation  regions  (13:19-20) 
were  found  to  correlate  with  two  rapid  expansions  of  the 
orbit  window. 
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Appendix 


c 

c 

c 

c 

c 

c 

c 

c 


Program  Phobos  written  by  Bob  Teets  Summer  1988 

integrates  a  version  of  Werner's  equations  of  motion  for 
a  massless  satellite  in  orbit  about  Phobos.  The  equations 
are  modified  by  the  introduction  of  the  system  of  units 


10 


15 


20 


25 


body  problem.  The  program  allows 
vector  to  data  files  for  plotting 
which  can  captured  in  a  log  file 


customary  for  the  three 
for  output  of  the  state 
or  simply  screen  output 
if  desired. 

common  /ham/  t , s ( 6 , 4 ) , ds ( 6 , 4 ) , err ( 6 ) , n , h 
double  precision  rO ( 6 ) , r ( 6 ) , dr ( 6 ) , drm ( 6 ) 
double  precision  t , s , ds , err , h , xtu 
logical  closed, plot , plane 
character*l  ans 
type  10 

f ormat ( 5x , ' Do  you  want  orbit  plots  ?  ',$ 


accept  1 5 , ans 
format ( al ) 

i f ( ( ans . eq . ' y ' ) . or . ( ans . eq , 


plot= . true . 
else 

plot= . false 
endif 
type  20 
format ( 5x, 'Do 
accept  15, ans 
if ( ( ans . eq . ' y ' ) 


If  so,  the  state  vector 
will  be  written  to  a 
file  state. dat;n,  where 
) )  then 

version,  n,  corresponds 
to  the  record  #  read 
from  the  file,  IC.dat, 
the  initial  conditions. 


you  want  plane  section  plots  ? 
or . ( ans . eq . ' Y ' 


) )  then 

If  so,  the  state  vector 
will  be  written  to  file 
section . dat ; n  twice  each 
orbit  as  y  passes  zero. 


status= ' old ' ) 
i=l , 6 ) , xtu 
input  initial 
in  Km,  m/sec. 


pos , 
T.U, 


vel,  tf 


6 ) , xtu 


echo  back  input 


plane= . true . 
else 

plane= . false . 
endif 

open ( uni t= 10 ,file='ic.dat' 
read ( 10 , 25 , end= 100 )  (r0(i) 

i 

format ( 7f 10 . 5 )  ! 

type  * , 'inputs  ',(r0(i),i=l 
do  i=l , 3 

rO ( i ) =r0 ( i ) /9 . 378d+03  !  converts  Km  to  D.U. 

r0(i+3)=r0(i+3)/9.378d+06*4.383d+03  !  converts 

enddo  !  m/sec  to  D.U. /T.U. 

do  i=l , 6 

s(i,l)=r0(i)  !  puts  initial  pos  and  vel  in  s 

dr { i ) =0 .  !  initialize  delta  state  vector 

drm(i)=0.  !  initialize  max  delta  state 

enddo 

s ( 4 , 1 ) =r0 ( 4 ) +1 . 0d+00-r0 ( 2 )  !  changes  xdot  to  Px 

s ( 5 , 1 ) =r0 ( 5 ) +r0 ( 1 )  !  changes  ydot  to  Py 


47 


norb=0 

nxt=0 

n=6 

h= 1 . 0d-04 
t=0 , Od+OO 
tlast=0 . 
tstep= . 01 
nsteps=xtu/h 
tp=0 . 

call  haming(nxt) 


!  orbit  counter 
!  haming  initialization  flag 
!  #  of  equations  of  motion 
!  delta  t,  time  step  for  integration 
!  set  initial  time  =  0 
!  time  at  end  of  last  orbit 
!  time  step  for  print  interval 
!  lets  integration  go  for  x  T.U. 

!  time  for  next  print 
!  initialize  haming 


if  (nxt  . eq .  0)  then 

type  *,  'haming  failed  to  initialize' 
stop 
endif 

if (plot )  open ( uni t=l 1 , f ile= ' state . dat ' , status= ' new ' ) 
if (plane) 

*  open (uni t= 12 , f ile= ' section . dat ' , status= ' new ' ) 
if(plot)  write(ll,30) 

*  t , (s ( i , nxt ) , i=l , 3) , (ds( i ,nxt ) , i=l , 3) 
format  (lx,7ell.4) 

do  it=l,nsteps 
do  i=l , 3 

r ( i ) =s ( i , nxt )  !  save  old  position  vector 

r ( i+3 ) =ds ( i f nxt )  !  save  old  velocity  vector 

enddo 

call  haming(nxt)  !  integrate  for  new 


!  integrate  for  new 
1  position  and  velocity 
!  print  state  vector 
!  set  next  print  time 


if(t.ge.tp)  then  !  print  state  vector 

tp=tp+tstep  !  set  next  print  time 

if(plot)  write(ll,30) 

*  t , (s( i ,  nxt ) , i=l , 3 ) , (ds ( i ,nxt ) , i=l , 3 ) 
endif 

r2=s ( 1 , nxt ) *s ( 1 , nxt )+s ( 2 , nxt ) *s ( 2 , nxt ) 

*  +s ( 3 , nxt ) *s ( 3 , nxt } 
if(r2.ge.l)  then 

type  orbit ', norb, '  beyond  1  D.U. ' 
if(plot)  write(ll,30) 

*  t , (s( i , nxt ) , i=l , 3) , (ds( i ,nxt ) , i=l , 3) 
if(plot)  close(ll) 

if(plane)  close(12) 

go  to  1  !  new  initial  conditions 

endif 

phobos=s ( 1 , nxt ) *s ( 1 , nxt ) /I . 326e-6 

*  +s ( 2 , nxt ) *s ( 2 , nxt ) /2 . 072e-6 

*  +s ( 3 , nxt ) *s ( 3 , nxt ) /I . 005e-6 
if (phobos . le . 1 . )  then 

type  orbit ', norb , '  ends  in  collision  with 

♦Phobos ' 

if(plot)  write(ll,30) 

*  t , ( s( i , nxt ) , i=l , 3) , (ds( i  ,nxt ) , i=l , 3 ) 
if(plot)  close(ll) 

if(plane)  close(12) 

go  to  1  !  new  initial  conditions 


j 


* 


* 


* 


* 


* 


* 

* 

* 

* 

* 


* 

* 

* 

* 

* 

* 


* 


endif 
do  i=l,3 

dr ( i ) =abs ( s ( i , nxt ) -r ( i ) )  !  change  in  position 

dr ( i+3 ) =abs (ds ( i , nxt ) -r ( i+3 ) )  !  change  in  velocity 

if (dr ( i ) . gt . drm( i ) )  drm(i)=dr(i)  !  max  delta  r 
if ( dr ( i+3 ) . gt . drm ( i+3 ) )  drm ( i+3 ) =dr ( i+3 )  !  max 

Idelta  v 

enddo 

if ( ( ( s ( 1 , nxt ) . It . 0 . 0d+00 ) . and. ( r ( 1 ) . gt . 0 . 0d+00) ) . or . 

( ( s ( 1 , nxt ) . gt . 0 . 0d+00 ) .and. (r(l) .lt.O.Od+OO) ) ) 
then 

type  *, 'passing  x  =  0,  xdot  =  1 , ds ( 1 , nxt ) , '  y  =  ' 

,  s  (  2  , nxt ) 

if (ds ( 1 , nxt ) . gt . 0 )  then 

norb=norb+l  !  increment  orbit  count 

period= (t-tlast)*4383. 
t last=t 

type  period  =', period,'  seconds  for  orbit' 

,  norb 
endif 

if(plot)  write(ll,30) 

t , ( s ( i , nxt ) ,i=l,3) , (ds(i, nxt ) , i=l , 3 ) 

endif 

if ( ( ( s { 2 , nxt ) . It . 0 . 0d+00 ) . and . ( r ( 2 ) . gt . 0 . 0d+00 ) ) . or . 

( ( s ( 2 , nxt ) .gt.0.0d+00) .and. (r(2) .lt.0.0d+00) ) ) 
then 

type  passing  y  =  0,  ydot  =  ' , ds ( 2 , nxt ) , '  x  =  ' 

if  (piot"Xtwrite  (11,30) 

t , ( s( i , nxt ) , i=l , 3 ) , (ds(i, nxt ) , i=l , 3 ) 
if(plane)  write(12,30) 

t , ( s ( i , nxt ) ,1=1,3) , (ds(i, nxt ) , i=l , 3 ) 

endif 

if ( ( ( s ( 3 , nxt ) . It . 0 . 0d+00 ) .and. (r(3) .gt.0.0d+00) ) .or. 

( ( s ( 3 , nxt ) . gt . 0 . 0d+00 ) .and. (r(3) . lt.0.0d+00) ) ) 
then 

type  passing  z  =  0,  zdot  =  ' , ds ( 3 , nxt ) , '  x  =  ' 
,s( 1 , nxt ) 

if(plot)  write(ll,30) 

t , ( s ( i , nxt ) , i=l , 3 ) , (ds(i, nxt ) ,1=1,3) 

endif 

if (norb . ge . 1 )  then  !  check  if  orbit  closed 

closed= . true . 
do  1=1,3 

if ( abs ( rO ( i ) -s { i , nxt ) ) . gt . drm( i ) ) 
closed= . false . 

if (abs ( rO ( i+3 ) -ds ( i , nxt ) ) . gt . drm ( i+3 ) ) 
closed= . false . 

enddo 

if (closed)  then 

type  *,'orbit  closed,  #  orbits  =  ',norb 
if(plot)  wrlte(ll,30) 

t , ( s ( i , nxt ) ,1=1,3) , (ds(i, nxt ) ,1=1,3) 
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!  new  initial  conditions 


100 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


if(plot)  close(ll) 
if(plane)  close{12 
go  to  1 
end  if 
endif 
enddo 

if(plot)  write(ll,30) 

t , { s ( i , nxt ) 
if(plot)  close{ll) 
if(plane)  close{12) 
go  to  1 
stop 
end 


) 


,1=1,3) , ( ds ( i , nxt ) , i=l , 3 ) 


!  new  initial  conditions 


subroutine  haming (nxt) 

haming  is  an  ordinary  differential  eqns  integrator 
it  is  a  fourth  order  predictor-corrector  algorithm 
which  means  that  it  carries  along  the  last  four 
values  of  the  state  vector,  and  extrapolates  these 
values  to  obtain  the  next  value  (the  prediction  part) 
and  then  corrects  the  extrapolated  value  to  find  a 
new  value  for  the  state  vector. 

the  value  nxt  in  the  call  specifies  which  of  the  4 
values  of  the  state  vector  is  the  "next"  one. 
nxt  is  updated  by  haming  automatically,  and  is  zero  on 
the  first  call 

the  user  supplies  an  external  routine  rhs(nxt)  which 
evaluates  the  equations  of  motion 

common  /ham/  x , s { 6 , 4 ) , ds ( 6 , 4 ) , errest ( 6 ) , n , h 
double  precision  x , s , ds , errest , h , hh , xo , tol 

all  of  the  good  stuff  is  in  this  common  block, 
x  is  the  independent  variable  (  time  ) 
s(6,4)  is  the  state  vector-  4  copies  of  it,  with  nxt 
pointing  at  the  next  one 

ds(6,4)  are  the  equations  of  motion,  again  four  copies 
a  call  to  rhs(nxt)  updates  an  entry  in  ds 
errest  is  an  estimate  of  the  truncation  error  - 
normally  not  used 

n  is  the  number  of  equations  being  integrated  -  6  here 
h  is  the  time  step 

tol  =  0 . 0000000001d+0 

switch  on  starting  algorithm  or  normal  propagation 
if (nxt )  190,10,200 

this  is  hamings  starting  algorithm ....  a  predictor  - 
corrector  needs  4  values  of  the  state  vector,  and  you 
only  have,  one-  the  initial  conditions. 
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c  haming  uses  a  Picard  iteration  (slow  and  painful)  to 

c  get  the  other  three. 

c  if  it  fails,  nxt  will  still  be  zero  upon  exit, 

c  otherwise  nxt  will  be  1,  and  you  are  all  set  to  go 

c 

10  xo  =  x 

hh  =  h/2 . 0d+00 
call  rhs ( 1 ) 
do  40  1  =  2,4 
x  =  x  +  hh 
do  20  i  =  1 , n 

20  s(i,l)  =  s(i,l-l)  +  hh*ds(i,l-l) 
call  rhs ( 1 ) 
x  =  x  +  hh 
do  30  i  =  1 ,  n 

30  s(i,l)  =  s(i,l-l)  +  h*ds(i,l) 

40  call  rhs ( 1 ) 
jsw  =  -10 
50  isw  =  1 

do  120  i  =  l,n 

hh  =  s ( i ,  1 )  +  h*  (  9 . 0d+00*ds ( i  ,  1  )  +  19 . 0d+00*ds ( i , 2  ) 
1  -  5 . 0d+00*ds ( i , 3 )  +  ds(i,4)  )  /  24.0d+00 

if(  dabs(  hh  -  s ( i , 2 ) )  .It.  tol  )  go  to  70 
isw  =  0 

70  s ( i , 2 )  =  hh 

hh  =  s ( i  ,  1 )  +  h*  (  ds ( i , 1 )  +  4 . 0d+00*ds ( i , 2 )  + 

*  ds ( i , 3 ) ) / 3  - 0d+00 

if(  dabs(  hh-s(i,3))  .It.  tol  )  go  to  90 
isw  =  0 

90  s ( i , 3 )  =  hh  hh=s ( i , 1 ) +h* ( 3 . 0d+00*ds ( i , 1 ) 

*  +9 . 0d+00*ds ( i , 2 ) +9 . 0d+00*ds (1,3) 

1  +  3 . 0d+00*ds (1,4)  )  /  8 . 0d+00 

if(  dabs ( hh-s ( i , 4 ) )  .It.  tol  )  go  to  110 
isw  =  0 

110  s(i,4)  =  hh 
120  continue 
x  =  xo 

do  130  1  =  2,4 
x  =  x  +  h 
130  call  rhs(l) 

if ( isw)  140,140,150 
140  jsw  =  jsw  +  1 

if(jsw)  50,280,280 
150  x  =  xo 
isw  =  1 
jsw  =  1 

do  160  i  =  1 ,  n 
160  errest(i)  =  0.0 
nxt  =  1 
go  to  280 
190  jsw  =  2 

nxt  =  labs (nxt) 
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c  this  is  hamings  normal  propagation  loop  - 

c 

200  x  =  X  +  h 

npl  =  mod{nxt,4)  +  1 
go  to  ( 2 10 , 230 ) , isw 
c  permute  the  index  nxt  modulo  4 

210  go  to  ( 270, 270 , 270 , 220 ), nxt 
220  isw  =  2 

230  nm2  =  mod(npl,4)  +  1 
nml  =  mod ( nm2 , 4 )  +  1 

npo  =  mod (nml, 4)  +  1 
c 

c  this  is  the  predictor  part 

c 

do  240  i  =  1 , n 

ds ( i , nm2 )  =  s(i.npl)  +  4  0d+00*h' I  2 . 0d+00*ds ( i , npo )  - 

1  ds(i.nml)  ♦  2 . Od+OO'ds ( i , nm? )  )  3 . 0d+00 

240  s(i,npl)  =  ds(i,nm2)  0 . 9256 1 9835*errest (  i  ) 
c 

c  now  the  corrector  -  fix  up  the  extrapolated  state 

c  based  on  the  better  value  of  the  equations  of  motion 

c 

call  rhs ( npl ) 
do  250  i  =  1  ,  n 

s(i,npl)  =  (  9 . 0d+00*s ( i , npo )  -  s ( i , nm2 )  +  3.0d+00*h*( 
1  ds(i,npl)  +  2 . 0d+00*ds ( i , npo )  -  ds(i,nml)))  /  8 . 0d+00 
errest(i)  =  ds(i,nm2)  -  s(i,npl) 

250  s ( i , npl )  =  s(i.npl)  +  0.0743801653  *  errest(i) 
go  to  ( 260 , 270 ) , jsw 
260  call  rhs(npl) 

270  nxt  =  npl 
280  return 
end 


subroutine  rhs (nxt) 

c 

c  rhs  contains  the  differential  equations  of  motion, 

c  the  basic  function  of  rhs  is  to  calculate  the 

c  equations  of  motion  (ds  =  f (s,t)  )  from  the  given 

c  current  state  (stored  in  s)  and  the  time  t. 

c  the  state  vector, s,  is  defined  as  follows: 

c  s( 1-3, nxt)  are  the  x,y,z  components  of  the  position 

c  vector, r.  s(4-6,nxt)  are  the  generalized  momenta,  Px, 

c  Py ,  Pz . 

c 

c  the  haming  common 

c 

common  /ham/  t , s ( 6 , 4 ) , ds ( 6 , 4 ) , err ( 6 ) , n , h 
double  precision  t,s,ds,err,h 
double  precision 

*  r,r2,r3,r4,r5,r7,x,x2,y,y2,z,z2,Px,Py,Pz 

double  precision 

r2x, r2y , r2z , xy , yz , Ixx , Iyy , Izz , I 1 , 12 , 13 , It , mO , ml 
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****************************************************** 


EVALUATE  THE  EQUATIONS  OF  MOTION 
****************************************************** 

ml  =  l .541036375d-08  !  M.U. 

rr.C  =  l  .  0d+00-ml  !  M.U. 

Ixx=ml*6. 153941736d-07  !  M.U.*D.U.2 

Iyy=ml *4 . 661904792d-07  !  M.U.*D.U.2 

Izz-ml *6 . 797057187d-07  !  M.U . *D .U. 2 

I 1=-Ixx+Iyy+Izz 

I2=Ixx-Iyy+Izz 

I3=Ixx+Iyy-Izz 

1 1 - 1 xx + 1 yy+ 1 z  z 

r2=s(l,nxt)*s(l,nxt)+s(2,nxt)*s(2,nxt) 

+s(3,nxt) *s ( 3 , nxt ) 
r=dsqrt ( r 2 ) 
r 3=r *r2 
r4  =  r2  *  r 2 
r5=r3*r2 
r7=r3*r4 

xy=s ( 1 , nxt ) *s ( 2 , nxt ) 
yz=s ( 2 , nxt ) *s ( 3 , nxt ) 
r2x=r2*s ( 1 , nxt ) 
r2y=r2*s ( 2 , nxt ) 
r2z-r2*s ( 3 , nxt ) 
x2=s ( 1 , nxt ) *s ( 1 , nxt ) 
y2=s ( 2 , nxt ) *s ( 2 , nxt ) 
z2=s(3,nxt) *s(3,nxt) 

ds(l,nxt)  =  s ( 4 , nxt ) +s ( 2 , nxt ) - 1 . Cd+00 
ds(2,nxt)  =  s { 5 , nxt ) -s ( 1 , nxt ) 
ds(3,nxt)  =  s(6,nxt) 

ds(4,nxt)  =  s ( 5 , nxt ) -mO* ( s ( 1 , nxt ) -3 . 0d+00*xy 

+1.5d+00*r2x)-ml*s(l,nxt)/r3+7.5d-01*It*s(l,nxt)/r5 
-7 . 5d-01 *s ( 1 , nxt ) /r7* ( ( 5 . 0d+00*x2-2 . 0d+00*r2 ) *11 
+5 . 0d+00*y2* i2+5 . 0d+00*z2*I3 ) 
ds(5.nxt)  =  -s ( 4 , nxt ) +m0* ( 1 . Od+OO-4 . 0d+00*s ( 2 , nxt ) 

-1 . 5d+00* ( r 2y-r2-2 . 0d+00*y2 ) ) -ml*s ( 2 , nxt ] /r3 
+  7 . 5d-0 1  *  1 1  *  s ( 2 , nxt) /r 5-7 . 5d-01*s ( 2 , nxt ) /r7* 

( ( 5 . 0d+00*y2-2 .*r2)*I2+5. 0d+00*x2*ll 
+5 . 0d+00*z2*I3 ) 

ds ( 6 , nxt )  =  -m0*(s(3,nxt)-3.0d+00*yz+1.5d+00*r2z)- 
ml*s(3,nxt)/r3+7.5d-01*It*s(3,nxt)/r5 
-7 . 5d-01*s ( 3 , nxt ) /r7* ( ( 5 . 0d+00*z2-2 . 0d+00*r2 ) *13 
+  5 . 0d+00 *x2  *11  +  5 . 0d+00*y2*I2 ) 
return 
end 
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